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We describe a formalism and numerical approach for studying spherically symmetric scalar field 
collapse for arbitrary spacetime dimension d and cosmological constant A. The presciption uses a 
double null formalism, and is based on field redefinitions first used to simplify the field equations in 
generic 2— dimensional dilaton gravity. The formalism is used to construct code in which d and A 
are input parameters. We reproduce known results in d = 4 and d = 6 with A = 0, and present new 



' results for d = 5 with zero and negative A. 
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^ ■ I. INTRODUCTION 

... ... .... . n . . 

^ ' It is an interesting fact that spherically symmetric gravitational collapse exhibits critical behaviour 111. This is a 

^^O I classical effect in phase space: there are one parameter families of initial data sets, for a variety of matter fields, such 
that as the parameter is tuned, a transition from reflection of infalling matter to black hole formation is observed 
numerically. Two types of behaviour are observed close to this transition point. Depending on the matter type, black 
holes form with zero or non-zero initial mass, and the matter field exhibits discrete or continuous self-similarity. A 

■ variety of matter fields have been studied since the initial seminal work by Choptuik. There have also been some 
, extensions beyond spherical symmetry, as well as a semi-analytic perturbation theory understanding of the critical 

■ exponent. (Recent reviews may be found in Refs. If an exact time dependent solution were available with the 
^H, appropriate boundary conditions, this critical behaviour would be manifested in the solution of the apparent horizon 
^ • condition, which is a transcendental equation [Q. However, to date only one time dependent (and non-self similar 

solution) is known, but this docs not have the required asymptotically flat boundary conditions Q. Thus the full 
• • I PDE problem must be studied numerically. 
. ^ ' Most numerical studies of spherically symmetric collapse have been in four spacetime dimensions. The only excep- 
: tions are the massless minimally coupled scalar field by Garfinkle Q in six spacetime dimensions, and by Pretorius 
^ and Choptuik and Husain and Olivier |^ in three dimensions with negative cosmological constant. In addition, 
5^ , there exist two papers that study this problem in the much simpler case of the self-similarity ansatz, an ODE problem, 
for any spacetime dimension |]lo| , |lT| . 

The purpose of this work is to present a formalism and numerical method for studying the gravitational collapse of 
a spherically symmetric scalar field for any value of spacetime dimension and cosmological constant. The approach 
reproduces and supplements known results in four and six dimensions with zero cosmological constant, and gives new 
results in five spacetime dimensions with zero and negative cosmological constant. 



II. FIELD EQUATIONS 

The basic idea for obtaining field equations valid for any dimension is to reduce the d— dimensional Einstein-scalar 
field equations by imposing spherical symmetry, and then use a field redefinition originally motivated by generic 2- 
dimensional dilaton gravity p^ . This allows the Einstein equations for d— dimensional, spherically symmetric scalar 
field collapse to be put into a form that can be managed numerically by a single code. In actuality, the formalism and 
code we describe below is applicable to a more general class of models that includes non-minimal scalar field coupling. 
This class of theories is quite broad, and contains as a sub-class the d— dimensional spherically symmetric case. 

Einstein gravity with cosmological constant in d spacetime dimensions is given by the action 
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The corresponding action for a minimally coupled scalar field is 

To impose spherical symmetry, we write the dimensional metric g^jy as 

dsld) = gapdx°'dx'^ + r'^{x°')dn(^a-2), 

where drJ(d-2) is the metric on S'''"^ and a, /3 = 1, 2. This gives the reduced action |] 

Stot '■— Sg + Sm 

= V(") / d'x^ ( ^^^^-^\R[g) - A) + \m + 



2G \^{n~l) 
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where n — d ~ 2^ V*-"-* is the volume of the unit n-sphere, and 

1 _ 8(n - 1) 
2G ~ IQim 
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The key simplification in our formalism is achieved by the following conformal reparametrization of the metric, 
which eliminates the kinetic term for (f) from the action |12| . Let 



where 



f)^((/.) = Cexp - 
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and G is an arbitrary constant. Now define a dimensionless "dilaton" field 
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Note that (j) is proportional to the area of the n-sphere at radius r. With these redefintions the reduced action takes 
the simpler form 



Stot = ^ / d^x^f^ 



0i?(.g) + y(")(0)l - / d'x^ H^''\<i))\dx? 



(12) 



where 



^In most cases one should perform dimensional reduction at the level of the field equations to guarantee that one obtains 
the correct solution space. It is well known |^] that in the present case the reduced field equations correspond to the field 
equations obtained from the reduced action. 
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and the overall factor of V^"^ has been dropped. 

For arbitrary functions V^(0) and H{if>), the action (|l2|) is that of generic dilat qn g ravity theory coupled to a scalar 
field in two spacetime dimensions. This theory has been studied in great detail [|l2| . The vacuum equations (x — 0) 
can be solved exactly. By choosing an adapted coordinate system in which the the dilaton (j) plays the role of the 
spatial coordinate (i.e. x = l(j)), the vacuum solution for the metric is 

ds^ = - 2GM)dt^ + {j{(t)) - 2GM)-^dx^ (15) 

where M plays the role of mass, and 

m= f dwC^) (16) 

Note that we have dropped the superscript n denoting spacetime dimension, since the above solution applies to the 
generic case. For the specific case of rf— dimensional spherically symmetric gravity 

n\ n J 

It can easily be verified by making the appropriate substitutions and conformal reparametrization (^ that the physical 
line element dsj^^-^ corresponding to ( p^ is precisely that of a d dimensional deSitter/anti-deSitter black hole with 
mass M.0 It is important to note that the metric (15) is singular at = even when M = 0. Up to numerical 



constants, j goes to zero as 

m ^ 01-1/" (18) 

near (f> — Q. This is not a physical singularity since the physical metric Ij is indeed the Minkowski metric when M and 
A are zero. Nonetheless, the vanishing of j(0) will affect the choice of boundary conditions in our numerical method. 

We now examine the field equations that derive from ( [T^ ) in double null coordinates, for which the metric may be 
parametrized as 

ds'^ — —2lg{u,v)(j)'{u,v)dudv (19) 

where the prime denotes partial differentiation with respect to the null coordinate v. (Recall that this is just the u — v 
part of the physical metric.) The corresponding field equations are 

0' = -^W")(0)#' (20) 

= 2G{xr (21) 



gif(»)(0) 

(F(")(0)x')- + (i?'"H'/')x)' = O. (22) 

In the above, the dot refers to differentiation with respect to u, which is treated like the "time" coordinate for the 
purposes of the following numerical integration. Remarkably, for arbitrary n Eqs. (|2C|-p2|) are virtually identical in 
form to those studied in ^ in the context of 2+1 dimensional AdS gravity.^ However, the boundary conditions are 



In order to get the overall scale factor right, one must choose the constant C in (^) appropriately. 
For n = 1 they are identical as expected. 
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special in the case of n = 1 (d = 3), so we will not consider this case further. Except where explicitly stated, we 
henceforth restrict consideration to n > 2 (which means spacetime dimension d > A). 

The evolution equations may be put in a form more useful for numerical solution by defining the variable 



h = X + 



(23) 



This effectively replaces the scalar field x by h. The evolution equations are 

= -5/2 



where 



and X is now to be considered a functional of h and g given by 



~g= / {g(^'V)dv'. 
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dv 



h(t>' 
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The integration constant K3{u) must be zero because the definition of h requires /i = x at = O.Q The function g is 
a functional of h and <j), obtained by integrating the constraint (|2l|): 



g = Ki (u) exp 



471 / dv^ih~x? 



(28) 



where Ki{u) is again an integration constant (i.e. independent of the "spatial" coordinate v). We consider the case 
of a spherically symmetric, collapsing shell of matter, with no black hole in the interior, initially. Thus, our boundary 
condition should be such as to guarantee that the metric at r = (which translates to ^ = 0) goes over to the vacuum 
solution (|l5|). By transforming the vacuum [M = 0) metric (|l^) to double null coordinates (|l^) we obtain a metric 
of the form 



where u = t — <b* and v = t - 



ds^ = —j{(j))dudv 
*, with the generalized "tortoise coordinate" (p* defined by 



(29) 
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d(f> 



With these definitions, 



/o 

corresponds to the surface v = u. Moreover, it follows that for the vacuum solution 
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Comparing metric ( p9[ ) to our general form (|T^) we see that g = 1 for the vacuum solution. Since we would like the 
numerical solution to approach the vacuum at = 0, the above analysis determines the required boundary conditions. 
In particular the integration constant Ki {u) = 1 , and 



which vanishes at (/> = in agreement with the expression (p6|). 



(32) 



*We will see below that although goes to zero at = it does so slowly enough to guarantee that the second term in ( [23 
vanishes at (/> = 0. 
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III. NUMERICAL METHOD 



The numerical scheme uses a v ('space') discretization to obtain a set of coupled ODEs: 

h{u,v) ^ hi{u), 4>{u,v) ^ (t>i{u). (33) 

where i = Q, ■ ■ ■ , N specifies the v grid. Initial data for these two functions is prescribed on a constant v slice, from 
which the functions g{u,v),g(u,v) are constructed. Evolution in the 'time' variable u is performed using the 4th. 
order Runge-Kutta method. The general scheme is similar to that used in [|l3|, together with some refinements used 
in ||]. This procedure was also used for the 3— dimensional collapse calculations in 

The initial scalar field configuration x(0, u = 0) is most conveniently specified as a function of (f> rather than r. 
(Recall that cj) cx r".) This together with the initial arrangement of the radial points 4>{v,u — 0) fixes all other 
functions. We used the initial specification 0(0, w) = v. 

We consider two types of initial scalar field configurations: the Gaussian and "tanh" functions 



Xg{u = 0, (/)) = at/) exp 

and 



(34) 



Xt{u = 0, 0) = a tanh(0). (35) 

These choices permit us to test "universality," which is defined to be the independence of the details of the collapse, 
such as the critical exponent, on the choice of initial data shapes and parameters. Although universality may be 
tested for a large variety of shapes and parameters, the emphasis here is on different spacetime dimensions, so we 
have restricted attention to varying the amplitude parameter for these two shapes. 

The initial values of the other functions are determined in terms of the above by computing the integrals for g„ 
and Qn using Simpson's rule. In all cases, we used values of 0o = 1 and cr = 0.3 for the Gaussian initial data. 

The boundary conditions at fixed u are 

0fe = O, 5fc = 0, gk = l. (36) 

where k is the index corresponding to the position of the origin (p — 0. (In the algorithm used, all grid points 
< i < k — 1 correspond to ingoing rays that have reached the origin and are dropped from the grid; see below). 
These conditions are equivalent to r{u,u) — 0, g\r=vi = 9{u^u) = 1, and guarantee regularity of the metric at r = 0. 
Notice that for our initial data, <j)k and hence hk are initially zero, and therefore remain zero at the origin because of 
Eqn. (H). 

As evolution proceeds via the Runge-Kutta procedure, the entries in the (pi array sequentially reach 0, at which 
point they are dropped from the grid. Thus the radial grid loses points with evolution. This is similar to the procedure 
used in |l| and |]. 

At each u step, a check is made to see if an apparent horizon has formed by observing the function 

ah = g^''d^<PdpP = -p, (37) 

whose vanishing signals the formation of an apparent horizon. ^ For each run of the code with fixed amplitude a, this 
function is scanned from larger to smaller radial values after each Runge-Kutta iteration, and evolution is terminated 
if the value of this function reaches 10~^. The corresponding radial coordinate value is recorded as Rah- In the 
subcritical case, it is expected that all the radial grid points reach zero without detection of an apparent horizon. 
This is the signal of pulse reflection. 

The results (a. Rah) are collated as in by seeking a relationship of the form 



^It is worth noting that there exist special foliations of black hole spacetimes which have no apparent horizon [|l^. Therefore 
there is in principle the possibility that a numerical scheme that encounters such a shcing may fail to detect black hole formation. 
This manifestly does not happen with the double null coordinates used in our study, as demonstrated by the form of the static 
solution described above, and our results below. 
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Rah oc (a - a*) 



(38) 



where a, is the critical amphtude which separates the black hole and reflection solutions. Rah is the basic dimensional 
scale if a black hole forms, and is linearly related to black hole mass in four spacetime dimensions. 

To improve numerical accuracy near (f) — Q we follow a procedure similar to 1^, where all functions on a constant 
u surface are expanded in power series in at = 0, and the first three values of the constraint integrals are derived 
using the respective power series. We write 

h = ha + hicj) (39) 

and calculate the parameters ft-Oi hi using the linear least squares fit for the first fifteen points in /i,i(u).^ From this 
the expressions for g and g follow. The remaining iV — 3 values of these functions are computed from their integral 
using Simpson's rule for equally spaced points. This linear fit near the </) = is necessary because it elegantly handles 
the problematic l/cj) factor in the h evolution equation, which would persist even if a finer mesh were used. The 
t)— derivatives of functions (needed for computing g and g) are calculated using f[ — {fi+i — fi-i)/2Av with end point 
values determined by linear extrapolation: /( = 2/2 — and = 2/^_;^ — /^_2- 

Further comments on the procedure are the following. The number of v grid points decreases as ingoing null 
geodesies cross r — 0, and so a refiected pulse cannot be followed back out toward infinity.^ Also, again due to the 
loss of grid points, and hence resolution, we are not able to observe the detailed behaviour of the scalar field very near 
criticality. An additional numerical adjustment concerns the enforcement of boundary conditions at the origin: there 
is the gravitational tendency for the matter to pile up at the origin as the collapse proceeds. However the formalism 
has the competing implicit condition h{u,u) = (ie. at = 0). This can lead to a shift of this boundary condition 
under evolution. It is rectified by adjusting the scalar field function at each time step by adding a constant shift at 
all points. This shift is of order 10~^ or less at each time step, and therefore there is a minor loss of accuracy, but a 
corresponding gain in stability. 

The code was tested for grid sizes ranging from 2000 to 6000 points, and with the u and v step sizes ranging from 
10^^ to 10~^, for the two types of initial data used, as well as the vacuum case of vanishing scalar field. These tests 
established that the code converges. Further tests of the code is the reproduction of the known results in four and six 
dimensions, which also demonstrates the accuracy of the results we obtain. All the results presented below were for 
a grid size of 6000 points, with u and v step sizes of 10""^. 

Finally, we point out that this procedure allows a more accurate analysis of the supercritical case than the subcritical 
one, again because of the loss of grid points as the evolution proceeds. This is less of an issue in the supercritical 
case because the number of points lost depends on the initial pulse amplitude, and our termination condition is such 
that typically one quarter of the original number of points are still present at termination, for the closest approach 
to criticality. We could of course monitor the subcritical case up to this point as well in an attempt to observe the 
self-similarity of the scalar field, but the very nature of this behaviour requires a way to replace lost grid points, as 
observed in |q]. For this reason we focus on supercritical evolutions. 



IV. RESULTS AND CONCLUSIONS 



The code was first tested to recalculate known results in four and six dimensions. It was then run for the five 
dimensional case with zero, positive and negative cosmological constant. All the calculations were performed for 
amplitudes above the threshold for black hole formation, and for initial data specified in both Gaussian and tanh 
forms, parametrized by amplitude a. The figures below show the scaling law Eqn. (|38|). The squares represent the 
points (o. Rah) and the lines are the least squares fit to these points. 

The 4— dimensional results for the Gaussian initial data are illustrated in Figure 1. The least squares fit gives a 
slope of 7 = 0.36, in good agreement (~ 4%) with earlier studies The figure also shows the oscillation about 

the fit line, again in accord with earlier work. 



There is nothing fundamental about this number, since the behaviour of the scalar field turns out to be very linear over the 
first several points; the results for ho and hi were virtually insensitive if the number of points used varied by a few on either 
side of fifteen. 

^A modification of our procedure along the lines suggested in may allow the tracking of the reflected pulse to future null 
infinity. 
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FIG. 1. Logarithmic plot of apparent horizon radius Rah versus initial scalar field amplitude {a — a*) in four spacetime 
dimensions for Gaussian initial data. 7 = 0.36 

The 6-dimensional results for Gaussian data appear in Figure 2. Our result for the critical exponent is 7 = 0.44. 
For comparison, the result in Rcf. (t) is 7 = 0.424. 
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FIG. 2. Logarithmic plot of apparent horizon radius Rah versus initial scalar field amplitude (a — a,) in six spacetime 
dimensions for Gaussian initial data. 7 = 0.44 

The results for the four and six dimensional calculations for the tanh initial data are 7 — 0.35 and 7 = 0.41 
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respectively. This provides a further check of our code, and further evidence of the insensitivity of the critical 
exponent to the shape of the initial data ("universality") in both four and six dimensions. Note that in six dimensions 
this is the first evidence for universality since Ref. contains results only for a specific gaussian form of initial data, 
different from the one used here. Note also that the agreement of our apparent horizon radius scaling results with 
those of the earlier works cited, shows that the apparent horizon appears to be a fairly good approximation to the 
event horizon of the long time static limit, (in sofar as these earlier studies actually find this limit). 

These tests of our formalism and code establish the consistency of our results with the earlier works mentioned 
above, and set the stage for new calculations for arbitrary values of the cosmological constant. Although our code 
allows calculations for any dimension, we focus on the 5— dimensional case mainly because results already exist for 3, 
4 and 6 dimensions. 

With A = 0, for the tanh initial data, we find a critical exponent of 7 = 0.41 (Figure 3). This value falls between 
the 4 and 6 dimensional cases as conjectured in W. 




With A = — 1, and Gaussian initial data, we find a critical exponent of 7 = 0.49 (Figure 4). All of the graphs 
show an oscillation about the least squares fit line. This is a known feature for zero cosmological constant, and is 
concomitant with discrete self-similarity of the critical solution. Our results for negative comological constant also 
show this feature, which indicates that the critical solution for this case also has discrete self-similarity. 

We find that in five dimensions the critical exponent appears not to be universal, at least in the supercritical 
approach to computing it. The Gaussian initial data yielded 7 — 0.52 for A = 0, in comparison to 7 = 0.41 for the 
tanh initial data (Figure 5). 



8 



i CO _ 



-10 



-8 -7 
ln(a-a*) 



FIG. 4. Logarithmic plot of apparent horizon radius Rah versus initial scalar field amplitude (a — a,) in five spacetime 
dimensions with A = — 1 for Gaussian initial data. 7 — 0.49 



ln(a-a*) 

FIG. 5. Logarithmic plot of apparent horizon radius Rah versus initial scalar field amplitude (a — a, 
dimensions with A = for Gaussian initial data. 7 — 0.52 



in five spacetime 



The reason for this is not clear to us and it would be worthwhile calculating the exponents using the subcritical 
approach suggested in Ref. |15|, where the Ricci scalar at the origin is calculated near criticality from below. It is also 
worth noting that a similar lack of universality is manifested in the 3— dimensional AdS case using the supercritical 
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apparent horizon method of computing 7 |16[ . 

We also carried out a preliminary investigation of the positive cosmological constant case in five dimensions. This 
is an interesting case because of the presence of a cosmological horizon in addition to the potential apparent horizon. 
Figure 6 shows graphs of the scalar field and apparent horizon functions h and ah in the left and right columns 
respectively, as functions of 0, prior to and at the onset of apparent horizon formation in the two successive rows. 
Note the location of the cosmological horizon in the right hand column near = 5.8. We find that the function 
v) evolved such that instead of the radial grid contracting as for the zero and negative A cases, it expanded as 
the scalar field moved towards the origin. This feature is visible in Figure 6: the range of in the lower graphes 
has expanded to 8 from 6. In fact the closer is the onset of apparent horizon formation, the larger the range of the 
(f> variable (and hence the radial grid). This prevented us from extracting accurate apparent horizon radii since the 
interesting features became confined to an ever shrinking part of the grid. We hope to study this in detail in future 
work. 





02468 02468 
phi phi 

FIG. 6. Plots of the scalar field h and the ah function Eqn. (^) prior to (top row) and nearer apparent horizon formation 
in five spacetime dimensions with positive A. Note the expansion of the (f> grid from 6 to 8 in the bottom two graphs. 
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In summary, we have described a formalism and code for studying spherically symmetric gravitational collapse of 
a scalar field for any d and A, presented new results in five dimensions, and given evidence for universality in six 
dimensions. In future work we will present results of a systematic analysis of the critical exponent as a function of 
both spacetime dimension and cosmological constant. 
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